home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / public / SciAn / src / ScianIsoCalcDelta.h < prev    next >
C/C++ Source or Header  |  1994-08-01  |  5KB  |  248 lines
  1. /*ScianIsoCalcDelta.h
  2.   Eric Pepke
  3.   23 March 1993
  4.  
  5.   Very evil include "macro" for calculating deltas in isosurfaces.  Delta is
  6.   not an accurate reflection of the gradient, but is the gradient normalized
  7.   to length one.  Since we're normalizing after the average, who cares?
  8.  
  9. Requires
  10.  
  11.     HEREC        coordinates here
  12.     HEREV        value here
  13.     FIELDHERE    field here
  14.     coord        coordinates field number
  15.     fieldNum    data field number
  16.  
  17. Optional
  18.     NEXTFIELD    field in +z
  19.     PREVFIELD    field in -z
  20.     PLUSIC        coordinates in plus i
  21.     MINUSIC        coordinates in minus i
  22.     PLUSJC        coordinates in plus j
  23.     MINUSJC        coordinates in minus j
  24.     PLUSKC        coordinates in plus k
  25.     MINUSKC        coordinates in minus k
  26. */
  27.  
  28. {
  29.     register long localOffset;
  30.     register real ds;        /*Change in scalar field*/
  31.     register real r;        /*Length of vector, squared*/
  32.     register real w;        /*Weight of vector*/
  33.     register real tx, ty, tz;    /*Total delta*/
  34.     register real dx, dy, dz;    /*Delta*/
  35.  
  36.     localOffset = index[0] * dims[1] + index[1];
  37.  
  38.     if (deltas[localOffset * 3] == missingData)
  39.     {
  40.     tx = 0.0;
  41.     ty = 0.0;
  42.     tz = 0.0;
  43.  
  44.     /*Minus i direction*/
  45.     if (index[0] > 0)
  46.     {
  47.         localOffset -= dims[1];
  48.         --index[0];
  49. #ifdef MINUSIC
  50.         dx = MINUSIC[0] - HEREC[0]; 
  51.         dy = MINUSIC[1] - HEREC[1]; 
  52.         dz = MINUSIC[2] - HEREC[2]; 
  53. #else
  54.         dx = SelectFieldComponent(coord, 0, index) - HEREC[0];
  55.         dy = SelectFieldComponent(coord, 1, index) - HEREC[1];
  56.         dz = SelectFieldComponent(coord, 2, index) - HEREC[2];
  57. #endif
  58.         ds = FIELDHERE[localOffset] - HEREV;
  59.  
  60.         ++index[0];
  61.         localOffset += dims[1];
  62.  
  63.         /*Weigh in component*/
  64.         r = (dx * dx + dy * dy + dz * dz);
  65.         if (r > 0.0)
  66.         {
  67.         w = ds / r;
  68.         tx += dx * w;
  69.         ty += dy * w;
  70.         tz += dz * w;
  71.         }
  72.     }
  73.  
  74.     /*Plus i direction*/
  75.     if (index[0] < dims[0] - 1)
  76.     {
  77.         localOffset += dims[1];
  78.         ++index[0];
  79. #ifdef PLUSIC
  80.         dx = PLUSIC[0] - HEREC[0]; 
  81.         dy = PLUSIC[1] - HEREC[1]; 
  82.         dz = PLUSIC[2] - HEREC[2]; 
  83. #else
  84.         dx = SelectFieldComponent(coord, 0, index) - HEREC[0];
  85.         dy = SelectFieldComponent(coord, 1, index) - HEREC[1];
  86.         dz = SelectFieldComponent(coord, 2, index) - HEREC[2];
  87. #endif
  88.         ds = FIELDHERE[localOffset] - HEREV;
  89.  
  90.         --index[0];
  91.         localOffset -= dims[1];
  92.  
  93.         /*Weigh in component*/
  94.         r = (dx * dx + dy * dy + dz * dz);
  95.         if (r > 0.0)
  96.         {
  97.         w = ds / r;
  98.         tx += dx * w;
  99.         ty += dy * w;
  100.         tz += dz * w;
  101.         }
  102.     }
  103.  
  104.     /*Minus j direction*/
  105.     if (index[1] > 0)
  106.     {
  107.         --localOffset;
  108.         --index[1];
  109. #ifdef MINUSJC
  110.         dx = MINUSJC[0] - HEREC[0]; 
  111.         dy = MINUSJC[1] - HEREC[1]; 
  112.         dz = MINUSJC[2] - HEREC[2]; 
  113. #else
  114.         dx = SelectFieldComponent(coord, 0, index) - HEREC[0];
  115.         dy = SelectFieldComponent(coord, 1, index) - HEREC[1];
  116.         dz = SelectFieldComponent(coord, 2, index) - HEREC[2];
  117. #endif
  118.         ds = FIELDHERE[localOffset] - HEREV;
  119.  
  120.         ++index[1];
  121.         ++localOffset;
  122.  
  123.         /*Weigh in component*/
  124.         r = (dx * dx + dy * dy + dz * dz);
  125.         if (r > 0.0)
  126.         {
  127.         w = ds / r;
  128.         tx += dx * w;
  129.         ty += dy * w;
  130.         tz += dz * w;
  131.         }
  132.     }
  133.  
  134.     /*Plus j direction*/
  135.     if (index[1] < dims[1] - 1)
  136.     {
  137.         ++localOffset;
  138.         ++index[1];
  139. #ifdef PLUSJC
  140.         dx = PLUSJC[0] - HEREC[0]; 
  141.         dy = PLUSJC[1] - HEREC[1]; 
  142.         dz = PLUSJC[2] - HEREC[2]; 
  143. #else
  144.         dx = SelectFieldComponent(coord, 0, index) - HEREC[0];
  145.         dy = SelectFieldComponent(coord, 1, index) - HEREC[1];
  146.         dz = SelectFieldComponent(coord, 2, index) - HEREC[2];
  147. #endif
  148.         ds = FIELDHERE[localOffset] - HEREV;
  149.  
  150.         --index[1];
  151.         --localOffset;
  152.  
  153.         /*Weigh in component*/
  154.         r = (dx * dx + dy * dy + dz * dz);
  155.         if (r > 0.0)
  156.         {
  157.         w = ds / r;
  158.         tx += dx * w;
  159.         ty += dy * w;
  160.         tz += dz * w;
  161.         }
  162.     }
  163.  
  164.     /*Minus k direction*/
  165.     if (index[2] > 0)
  166.     {
  167.         --index[2];
  168.  
  169. #ifdef MINUSKC
  170.         dx = MINUSKC[0] - HEREC[0];
  171.         dy = MINUSKC[1] - HEREC[1];
  172.         dz = MINUSKC[2] - HEREC[2];
  173. #else
  174.         dx = SelectFieldComponent(coord, 0, index) - HEREC[0];
  175.         dy = SelectFieldComponent(coord, 1, index) - HEREC[1];
  176.         dz = SelectFieldComponent(coord, 2, index) - HEREC[2];
  177. #endif
  178.  
  179.  
  180. #ifdef PREVFIELD
  181.         ds = PREVFIELD[localOffset] - HEREV;
  182. #else
  183.         ds = SelectFieldScalar(fieldNum, index) - HEREV;
  184. #endif
  185.         ++index[2];
  186.  
  187.         /*Weigh in component*/
  188.         r = (dx * dx + dy * dy + dz * dz);
  189.         if (r > 0.0)
  190.         {
  191.         w = ds / r;
  192.         tx += dx * w;
  193.         ty += dy * w;
  194.         tz += dz * w;
  195.         }
  196.     }
  197.  
  198.     /*Plus k direction*/
  199.     if (index[2] < dims[2] - 1)
  200.     {
  201.         ++index[2];
  202. #ifdef PLUSKC
  203.         dx = PLUSKC[0] - HEREC[0]; 
  204.         dy = PLUSKC[1] - HEREC[1]; 
  205.         dz = PLUSKC[2] - HEREC[2]; 
  206. #else
  207.         dx = SelectFieldComponent(coord, 0, index) - HEREC[0];
  208.         dy = SelectFieldComponent(coord, 1, index) - HEREC[1];
  209.         dz = SelectFieldComponent(coord, 2, index) - HEREC[2];
  210. #endif
  211.  
  212. #ifdef NEXTFIELD
  213.         ds = NEXTFIELD[localOffset] - HEREV;
  214. #else
  215.         ds = SelectFieldScalar(fieldNum, index) - HEREV;
  216. #endif
  217.  
  218.         --index[2];
  219.  
  220.         /*Weigh in component*/
  221.         r = (dx * dx + dy * dy + dz * dz);
  222.         if (r > 0.0)
  223.         {
  224.         w = ds / r;
  225.         tx += dx * w;
  226.         ty += dy * w;
  227.         tz += dz * w;
  228.         }
  229.     }
  230.  
  231.     /*Put it all together*/
  232.     {
  233.         real r;
  234.         r = sqrt(tx * tx + ty * ty + tz * tz);
  235.         if (r > 0.0)
  236.         {
  237.         r = 1.0 / r;
  238.         tx *= r;
  239.         ty *= r;
  240.         tz *= r;
  241.         }
  242.         deltas[localOffset * 3] = tx;
  243.         deltas[localOffset * 3 + 1] = ty;
  244.         deltas[localOffset * 3 + 2] = tz;
  245.     }
  246.     }
  247. }
  248.